% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme choses the preference parameters (CRRA rho,
% discount factor beta) according to various criteria, and produces tables
% and graphs

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version September 2021
% =====================================================================%
clear
clf
% this should be the replications folder
Tessa1Tobi2=2; % set to 1 for tessa, 2 for tobi
if Tessa1Tobi2==1
    currentfolder='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication';
elseif Tessa1Tobi2==2
    currentfolder='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication';
end
COALDEV=[1,0,2]%,2,3,4]; % benchmark: [1,0,2] %0 for ID, 1 for CD, 2 for SI, 3 for full insurance, 4 for heterogeneous preferences ID
Rov_sb=0; % benchmark=0; outside option for ROV, benchmark: 0
one_per_aut=0; % benchmark=0; one-period autarky, benchmark: 0
fixedeffect=1; % benchmark=1; fixed effects in model and data, benchmark: 1
unconditional=0; % benchmark: 0 % set to 1 if you want to use raw data for cons-inc moments, and model output with income process based on raw data
rhogrid=0; % benchmark: 0 % set this to 1 if you want to draw a picture for combinations of rho, beta
est_rho1=0;  % benchmark: 0   % estimate the AR(1) coefficient, rather than the measurement error
Nincproc=10; % benchmark=1% set to 10 for meas error; number of income processes
ext=1; % put to 1 if you want to make a table with CD, ID and SI
pres=0;   % set to 1 for presentation tables
phicalib=0; % set to nan if you want to estimate the borrowing limit, otherwise set to 0
if pres==0
    tablefontsize='small';
else
    tablefontsize='tiny';
end
buildup=0; % set to 1 if you're dealing with buildup files
buildupvec=[0,4];
samesize=0; % set to 1 if you want to consider the ID economy of the same size as the CD economy
samedelta=0;
deltamin=0.5;
rho1min=-0.6;
% input and output directories etc
write=1 ;% set to 1 if you want to write tables
analytic=1; % set to 1 if you want to use the analytic variance for the regression coefficients (all other moments and covariances are bootstrapped)
outlier=0; % set to 0 for no outliers removed, 1 for largest observations in village 1 in year 2, 2 for all of round 2 removed
beta_calib=0.9; % a priori value for beta
rho_calib=1; %... and for CRRA
bootstrapSE=2; % set this to 1 for bootstrapped standard errors of the relative regression coefficients and variances; set to 2 for analytical / STATA standard errors
Date='_10_Sep_2021';
printtable=0;
relvarincscale=1;
Momvec1data=NaN(18,1);
Momvec2data=NaN(18,1);

% =====================================================================%
% ====    Data moments
% =====================================================================%

clear icrisatdata momentsdatabootstrap DCdata DYdata DCdataneg DYdatapos DCdataneg DYdataneg Cdata Ydata VCV1data VCV2data
clear RelCYVardata RelDCDYVardata betaDCDYdatadyneg betaDCDYdatadypos betaDCDYdata Momvec1datamat Momvec2datamat DCskewdata RelDCDYVardata Cdata Ydata DCdata DYdata
clear betaDCDYdata CVARdata YVARdata DCVARdata DYVARdata Cskewdata DCskewdata RelCYVardata RelDCDYVardata

if buildup==1
    VILLAGE=1;
elseif Rov_sb==1;
    VILLAGE=[1,2,3];
elseif rhogrid==1;
    VILLAGE=[1];
else
    VILLAGE=[1,2,3];
end

if Tessa1Tobi2==1
    outputpath='..\Paper\Latex';
    if unconditional==1 && one_per_aut==0
        Date=[Date '_uncond_no_one_per_aut'];
        modeldatapath=[currentfolder '\Programmes\Main Results and Robustness'];
    else
        modeldatapath=[currentfolder '\Programmes\Main Results and Robustness'];
    end
    
    if  est_rho1==1
        Date=[Date '_estrho1'];
    elseif rhogrid==1
        Date=[Date '_rhogrid'];
    elseif  Rov_sb==1
        Date=[Date '_Rov_sb'];
    elseif one_per_aut==1
        Date=[Date '_cond'];
    elseif fixedeffect==0
        Date=[Date '_no_one_per_aut_no_FE'];
    elseif unconditional==0 && one_per_aut==0
        Date=[Date '_no_one_per_aut'];
    end
    modeldatapathvil1=modeldatapath;
    datapath=[currentfolder '\Data\Output'];
    datapathvil1=datapath;
    addpath('.');
else
    outputpath='..\Paper\Latex';
    if unconditional==1 && one_per_aut==0
        Date=[Date '_uncond_no_one_per_aut'];
        modeldatapath=[currentfolder '\Programmes\Main Results and Robustness'];
    else
        modeldatapath=[currentfolder '\Programmes\Main Results and Robustness'];
    end
    
    if  est_rho1==1
        Date=[Date '_estrho1'];
    elseif rhogrid==1
        Date=[Date '_rhogrid'];
    elseif  Rov_sb==1
        Date=[Date '_Rov_sb'];
    elseif one_per_aut==1
        Date=[Date '_cond'];
    elseif fixedeffect==0
        Date=[Date '_no_one_per_aut_no_FE'];
    elseif unconditional==0 && one_per_aut==0
        Date=[Date '_no_one_per_aut'];
    end
    modeldatapathvil1=modeldatapath;%sprintf('C:/Users/tbroer/Dropbox/Research/Current_Projects/tessa/Programmes/Output estimation 15 Dec outlier %d', outlier);
    datapath=[currentfolder '\Data\Output'];
    datapathvil1=datapath;
    addpath('.');
end

for village=VILLAGE
    datafile=['icrisatmoments' int2str(outlier) ];
    
    if village==1
        cd(datapathvil1)
    else
        cd(datapath)
    end
    eval(['load ' datafile '.out'])
    eval(['datamoments=' datafile ';']);
    eval(['clear ' datafile])
    Nom=4; %no. of moments
    analyticmoments=datamoments(1:Nom,1:3);
    bootstrapvar=datamoments(1:Nom,4+(village-1)*Nom+1:4+(village)*Nom);
    analyticvar=datamoments(1:Nom,16+(village-1)*Nom+1:16+(village)*Nom);
    %%%in case you want to use the exact variance for the regression
    %%%coefficients
    analyticrawmoments=datamoments(1:Nom,28+1:28+3);
    bootstraprawvar=datamoments(1:Nom,28+4+(village-1)*Nom+1:28+4+(village)*Nom);
    analyticrawvar=datamoments(1:Nom,28+16+(village-1)*Nom+1:28+16+(village)*Nom);
    %%%in case you want to use the exact variance for the regression
    %%%coefficients
    
    VCV2conddata=bootstrapvar;
    VCV2rawdata=bootstraprawvar;
    
    if analytic==1
        VCV2conddata(2,2)=analyticvar(2,2);
        VCV2conddata(4,4)=analyticvar(4,4);
        VCV2rawdata(2,2)=analyticrawvar(2,2);
        VCV2rawdata(4,4)=analyticrawvar(4,4);
    end
    
    Momvec2conddata=analyticmoments(:,village)';
    Momvec2rawdata=analyticrawmoments(:,village)';
    if unconditional==1
        VCV2data=VCV2rawdata;
        Momvec2data=Momvec2rawdata;
    else
        VCV2data=VCV2conddata;
        Momvec2data=Momvec2conddata;
    end
    
    COUNT=0;
    Countadd=0;
    for coaldev=COALDEV
        COUNT=COUNT+1;
        if village==1
            cd(modeldatapathvil1)
            if buildup==0
                maxvss=33;
            elseif buildup==1
                maxvss= 33;
            end
        elseif village==2
            cd(modeldatapath)
            if buildup==0
                maxvss=36;
            elseif buildup==1
                maxvss=36;
            end
        else
            cd(modeldatapath)
            if buildup==0
                maxvss=30;
            elseif buildup==1
                maxvss=36;
            end
        end
        %loading in once to get the dimensions
        if est_rho1==1
            eval(['load EstrhoMomentsest' int2str(village) int2str(coaldev) '10' int2str(fixedeffect) int2str(one_per_aut) '.mat'])
        elseif buildup== 1 && coaldev==0
            eval(['load BuildupMomentsest' int2str(village) int2str(coaldev) '10' int2str(fixedeffect) int2str(one_per_aut) '.mat'])
        elseif Rov_sb==1
            clear Momentsest
            eval(['load SBMomentsest' int2str(village) int2str(coaldev) '10' int2str(fixedeffect) int2str(one_per_aut) '.mat'])
        elseif rhogrid==1
            clear Momentsest
            eval(['load GridMomentsest' int2str(village) int2str(coaldev) '10.mat'])
        elseif coaldev==4
            clear Momentsest
            eval(['load Momentsest' int2str(village) int2str(coaldev) '10.mat'])
        elseif one_per_aut==1
            clear Momentsest
            eval(['load Momentsest' int2str(village) int2str(coaldev) '10' int2str(fixedeffect) int2str(one_per_aut) '.mat'])
        else
            eval(['load No_1_per_aut_Momentsest' int2str(village) int2str(coaldev) '1' int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut) '.mat'])
        end
        dim1= size(Momentsest,2); dim2=size(Momentsest,3); dim3=size(Momentsest,4);
        clear W4 W2 W4d W2d W2m W4dm W2dm W4dmr W2dmr Momentstable Momentsrawtable Momentssum
        W4=ones(31,dim2,dim3,maxvss,Nincproc)*Inf; W2=W4; W4d=W4; W2d=W4;
        W4m=ones(31,dim2,dim3,1,Nincproc)*Inf; W2m=W4m; W4dm=ones(31,dim2,dim3,Nincproc)*Inf; W2dm=W4m;
        Momentstable=nan(10,dim1,dim2,dim3,maxvss,Nincproc);
        Momentsrawtable=nan(9,dim1,dim2,dim3,maxvss,1);
        Momentssum=[];
        
        clear Momentsest
        for incproc=1:Nincproc % loop over different income processes implied by income measurement error
            if est_rho1==1
                eval(['load EstrhoMomentsest' int2str(village) int2str(coaldev) int2str(incproc)  '0' int2str(fixedeffect) int2str(one_per_aut) '.mat'])
            elseif buildup==1 && coaldev==0
                eval(['load BuildupMomentsest' int2str(village) int2str(coaldev) '10' int2str(fixedeffect) int2str(one_per_aut) '.mat'])
            elseif rhogrid==1
                clear Momentsest
                eval(['load GridMomentsest' int2str(village) int2str(coaldev) '10.mat'])
            elseif Rov_sb==1
                clear Momentsest
                eval(['load SBMomentsest' int2str(village) int2str(coaldev) '10' int2str(fixedeffect) int2str(one_per_aut) '.mat'])
            elseif coaldev==4
                clear Momentsest
                eval(['load Momentsest' int2str(village) int2str(coaldev) int2str(fixedeffect) int2str(one_per_aut) '.mat'])
            elseif one_per_aut==1
                clear Momentsest
                eval(['load Momentsest' int2str(village) int2str(coaldev) '10' int2str(fixedeffect) int2str(one_per_aut) '.mat'])
            else
                eval(['load No_1_per_aut_Momentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut) '.mat'])
            end
            vssset=[find(isfinite(squeeze(Momentsest(1,1,end,1,:))) & squeeze(Momentsest(1,1,end,1,:))>0)];
            
            if isempty(vssset)==1
                vssset=[find(isfinite(squeeze(Momentsest(1,1,1,1,:))) & squeeze(Momentsest(1,1,1,1,:))>0)];
                rhoset=[find(isfinite(squeeze(Momentsest(1,1,1,:,vssset(1)))) & squeeze(Momentsest(1,1,1,:,vssset(1)))>0)]'  ;
            else
                rhoset=[find(isfinite(squeeze(Momentsest(1,1,end,:,vssset(1)))) & squeeze(Momentsest(1,1,end,:,vssset(1)))>0)]'  ;
                if coaldev==2 && isnan(phicalib)==0
                    rhoset=find(squeeze(Momentsest(29,1,end,:,vssset(1)))==phicalib);
                end
            end
            if buildup==1 && coaldev==0 && samesize==1
                vssset=vssset(vssset==(solution_4(1,COUNT-1)-1))+buildupvec;
            end
            
            %%%get out betavec and rhovec
            if  incproc==1
                betavec=[squeeze(Momentsest(2,1,:,1,vssset(1),1))];
                rhovec=[squeeze(Momentsest(3,1,1,:,vssset(1),1))];
            end
            
            % ================
            % vector of moments
            % ================
            % col 1-5 are parameters
            % 6-11 relative variance dc/dy, dc(dy>0)/dc(dy<0);  unconditional (6:7), conditional on village income
            % (8:9), conditonal on group income (10:11)
            % 12-14: beta, beta high, beta low unconditional
            % 15-17: beta, beta high, beta low conditional on
            % village income
            % 18-20: beta, beta high, beta low conditional on
            % group income
            % 21:22 are variances of dclog, dylog
            if coaldev==4
                deltaset=[find(isfinite(squeeze(Momentsest(1,1,:,1,vssset(1)))) & squeeze(Momentsest(1,1,:,1,vssset(1)))>0  & squeeze(Momentsest(2,1,:,1,vssset(1)))>=deltamin) ]';
            else
                deltaset=[find(isfinite(squeeze(Momentsest(1,1,:,1,vssset(1)))) & squeeze(Momentsest(1,1,:,1,vssset(1)))>0  & squeeze(Momentsest(2,1,:,1,vssset(1)))>=deltamin) ]';
            end
            deltaset
            for inddelta=deltaset
                for indrho=rhoset
                    if sum((squeeze(Momentsest(1,1,inddelta,indrho,:)))>0 & isfinite(squeeze(Momentsest(1,1,inddelta,indrho,:))))>0
                        for j=1:size(Momentsest,2)
                            vss_set=[find((squeeze(Momentsest(1,j,inddelta,indrho,:)))>0 & isfinite(squeeze(Momentsest(1,j,inddelta,indrho,:))))]';
                            vss_max=max(vss_set);
                            if coaldev==0
                                if buildup==0
                                    vss_set=vss_max;
                                elseif buildup==1 && samesize==1
                                    vss_set=solution_4(1,COUNT-1)-1+buildupvec;
                                end
                            end
                            for vss=vss_set
                                
                                Vardylog(village)=Momentsest(22,j,inddelta,indrho,vss);
                                Vardylog_cond(village)=Momentsest(30,j,inddelta,indrho,vss);
                                
                                Momentstable(1:3,j,inddelta,indrho,vss,incproc)=Momentsest(1:3,j,inddelta,indrho,vss);
                                Momentstable(4:7,j,inddelta,indrho,vss,incproc)=[Momentsest([8 15],j,inddelta,indrho,vss) ; Momentsest(9,j,inddelta,indrho,vss)/Vardylog_cond(village) ;(Momentsest([16],j,inddelta,indrho,vss)-Momentsest([17],j,inddelta,indrho,vss))];
                                Momentstable(8,j,inddelta,indrho,vss,incproc)=(Momentsest(21,j,inddelta,indrho,vss)-Momentsest(21,1,inddelta,indrho,vss))./Momentsest(21,j,inddelta,indrho,vss);%2*maxcerr*(j-1)/(size(moments,2)-1)/moments(28,j);
                                if est_rho1==0
                                    Momentstable(9,j,inddelta,indrho,vss,incproc)=2*Momentsest(5,j,inddelta,indrho,vss)./Vardylog(village);%Momentsest(22,j,inddelta,indrho,vss);
                                elseif est_rho1==1
                                    Momentstable(9,j,inddelta,indrho,vss,incproc)=Momentsest(4,j,inddelta,indrho,vss);
                                    if Momentstable(9,j,inddelta,indrho,vss,incproc)<rho1min
                                        Momentstable(4:8,j,inddelta,indrho,vss,incproc)=nan;
                                    end
                                end
                                %actual group size rather than size of rest
                                %of village
                                Momentstable(1,j,inddelta,indrho,vss,incproc)=Momentstable(1,j,inddelta,indrho,vss,incproc)+1;
                                
                                if coaldev==2
                                    % put this if you want to have the borrowing limit
                                    % Momentstable(10,j,inddelta,indrho,vss,incproc)=Momentsest(29,j,inddelta,indrho,vss);
                                end
                                
                                Momentsrawtable(1:3,j,inddelta,indrho,vss,incproc)=Momentsest(1:3,j,inddelta,indrho,vss);
                                Momentsrawtable(4:7,j,inddelta,indrho,vss,incproc)=[Momentsest([6 12],j,inddelta,indrho,vss) ;Momentsest(7,j,inddelta,indrho,vss)/Vardylog(village); Momentsest([13],j,inddelta,indrho,vss)-Momentsest([14],j,inddelta,indrho,vss)];
                                Momentsrawtable(8,j,inddelta,indrho,vss,incproc)=(Momentsest(21,j,inddelta,indrho,vss)-Momentsest(21,1,inddelta,indrho,vss))./Momentsest(21,j,inddelta,indrho,vss);%2*maxcerr*(j-1)/(size(moments,2)-1)/moments(28,j);
                                if est_rho1==0
                                    Momentsrawtable(9,j,inddelta,indrho,vss,incproc)=2*Momentsest(5,j,inddelta,indrho,vss)./Momentsest(22,j,inddelta,indrho,vss);%Momentsest(22,j,inddelta,indrho,vss);
                                elseif est_rho1==1
                                    Momentsrawtable(9,j,inddelta,indrho,vss,incproc)=Momentsest(4,j,inddelta,indrho,vss);
                                end
                                Momentsrawtable(1,j,inddelta,indrho,vss,incproc)=Momentsrawtable(1,j,inddelta,indrho,vss,incproc)+1;
                                
                                MOMENTSDIAGNOSTICS(1:9,inddelta,indrho,vss,COUNT,village)=Momentstable(1:9,1,inddelta,indrho,vss,1);
                                if j==1
                                    MOMENTSDIAGNOSTICS(10,inddelta,indrho,vss,COUNT,village)=Momentsest(10,1,inddelta,indrho,vss);
                                end
                                momentsdev=Momentstable(4:7,j,inddelta,indrho,vss,incproc)- Momvec2data(1,1:4)';
                                %%%2,4, d:diagonal, m: max stable group/village
                                %%% r=raw
                                
                                W4(j,inddelta,indrho,vss,incproc)=momentsdev'*inv(VCV2data(1:4,1:4))*momentsdev;
                                W4d(j,inddelta,indrho,vss,incproc)=momentsdev'*inv(diag(diag(VCV2data(1:4,1:4))))*momentsdev;
                                
                                W2(j,inddelta,indrho,vss,incproc)=momentsdev(1:2)'*inv(VCV2data(1:2,1:2))*momentsdev(1:2);
                                W2d(j,inddelta,indrho,vss,incproc)=momentsdev(1:2)'*inv(diag(diag(VCV2data(1:2,1:2))))*momentsdev(1:2);
                            end
                            %this keeps the result for the largest group
                            %size only
                            if isempty(vss_set)==0
                                W4m(j,inddelta,indrho,incproc)=W4(j,inddelta,indrho,vss,incproc);
                                W4dm(j,inddelta,indrho,incproc)=W4d(j,inddelta,indrho,vss,incproc);
                                W2m(j,inddelta,indrho,incproc)=W2(j,inddelta,indrho,vss,incproc);
                                W2dm(j,inddelta,indrho,incproc)=W2d(j,inddelta,indrho,vss,incproc);
                            end
                        end
                        if Momentsest(2,1,inddelta,indrho,vss)>beta_calib-0.000001 & Momentsest(2,1,inddelta,indrho,vss)<beta_calib+0.000001 & Momentsest(3,1,inddelta,indrho,vss)==rho_calib & incproc==1 && ( coaldev~=2 || Momentsest(29,1,inddelta,indrho,vss)==0)
                            Momentscalib(1:size(Momentsest,1),COUNT)=Momentsest(:,1,inddelta,indrho,vss)+[1;zeros(size(Momentsest,1)-1,1)];
                        end
                        
                        if incproc==1 & isempty(vss_set)==0
                            Momentssum=[Momentssum Momentstable(:,1,inddelta, indrho,vss,incproc)];
                        end
                    end
                end
            end
            incproc
        end
        
        % ================
        % Estimation
        % ================
        if rhogrid==1
            for Qbeta=1:length(deltaset)
                if coaldev==1
                    vss_max_plot(coaldev+1,Qbeta)=max(find(isfinite(squeeze(Momentsest(1,1,Qbeta,1,:))).*(~isnan(squeeze(Momentsest(1,1,Qbeta,1,:)))).*squeeze(Momentsest(1,1,Qbeta,1,:))>0))
                    vss_max_plot2m(coaldev+1,Qbeta)=vss_max_plot(coaldev+1,Qbeta);
                    max4(coaldev+1)=Qbeta;
                elseif coaldev==0
                    vss_max_plot(coaldev+1,Qbeta)=max(Momentsest(1,1,1,1,:));
                    vss_max_plot2m(coaldev+1,Qbeta)=vss_max_plot(coaldev+1,Qbeta);
                    max2(coaldev+1)=Qbeta;
                end
                CC(Qbeta)=(squeeze(W4dm(1,Qbeta,1,1))') ;
                CC2m(Qbeta)=(squeeze(W2dm(1,Qbeta,1,1))') ;
                beta_plot(Qbeta,coaldev+1,village)=Momentsest(2,1,Qbeta,1,vss_max_plot(coaldev+1,Qbeta))
                sigma_plot(Qbeta,coaldev+1,village)=Momentsest(3,1,Qbeta,1,vss_max_plot(coaldev+1,Qbeta))
                moment_plot(Qbeta,1,coaldev+1,village)=Momentsest(8,1,Qbeta,1,vss_max_plot(coaldev+1,Qbeta))
                moment_plot(Qbeta,2,coaldev+1,village)=Momentsest(15,1,Qbeta,1,vss_max_plot(coaldev+1,Qbeta))
                moment_plot(Qbeta,3,coaldev+1,village)=Momentsest(9,1,Qbeta,1,vss_max_plot(coaldev+1,Qbeta))/Vardylog_cond(village);
                moment_plot(Qbeta,4,coaldev+1,village)=Momentsest(16,1,Qbeta,1,vss_max_plot(coaldev+1,Qbeta))-Momentsest(17,1,Qbeta,1,vss_max_plot(coaldev+1,Qbeta))
                data_plot(Qbeta,:,coaldev+1,village)=Momvec2data;
                fit_plot(Qbeta,coaldev+1,village)=CC(Qbeta);
                Sigma_plot(village,Qbeta,coaldev+1,village)=sigma_plot(Qbeta,coaldev+1,village);
                
                sigma_plot2m(Qbeta,coaldev+1,village)=Momentsest(3,1,Qbeta,1,vss_max_plot2m(coaldev+1,Qbeta))
                moment_plot2m(Qbeta,1,coaldev+1,village)=Momentsest(8,1,Qbeta,1,vss_max_plot2m(coaldev+1,Qbeta))
                moment_plot2m(Qbeta,2,coaldev+1,village)=Momentsest(15,1,Qbeta,1,vss_max_plot2m(coaldev+1,Qbeta))
                moment_plot2m(Qbeta,3,coaldev+1,village)=Momentsest(9,1,Qbeta,1,vss_max_plot2m(coaldev+1,Qbeta))/Vardylog_cond(village);
                moment_plot2m(Qbeta,4,coaldev+1,village)=Momentsest(16,1,Qbeta,1,vss_max_plot2m(coaldev+1,Qbeta))-Momentsest(17,1,Qbeta,1,vss_max_plot2m(coaldev+1,Qbeta))
                data_plot2m(Qbeta,:,coaldev+1,village)=Momvec2data;
                fit_plot2m(Qbeta,coaldev+1,village)=CC2m(Qbeta);
                Sigma_plot2m(village,Qbeta,coaldev+1,village)=sigma_plot2m(Qbeta,coaldev+1,village);
            end
            if coaldev==COALDEV(end)
                try
                    identif_pic
                end
            end
        else
            if coaldev<2% only do if ID or CD, not with savings
                if coaldev==0 % initialise some variables
                    delta_4_gs_ms=nan(9,3,2);delta_4_gs=nan(9,1,2);delta_2_gs_ms=nan(9,3,2);delta_2_gs=nan(9,1,2);
                end
                % ===============
                %four moments
                % ===============
                if Nincproc>1
                    WW=W4d(:,deltaset,rhoset,vssset,:);
                    criterion=reshape(WW,prod(size(WW)),1);
                    [CC,II]=min(criterion);
                    [II1,II2,II3,II4,II5]=ind2sub(size(WW),II);
                    solution_4_gs_ms(:,COUNT)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5); CC];
                    if II2<length(deltaset) && II2>1
                        delta_4_gs_ms(:,1,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))];
                    elseif II2==length(deltaset)
                        delta_4_gs_ms(:,1,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))];
                    elseif II2==1
                        delta_4_gs_ms(:,1,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))];
                    end
                    if II1<size(W4d,1) && II1>1
                        delta_4_gs_ms(:,2,COUNT)=[(Momentstable(1:9,II1-1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1+1,deltaset(II2),rhoset(II3),vssset(II4),II5))/(Momentstable(8,II1-1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(8,II1+1,deltaset(II2),rhoset(II3),vssset(II4),II5))];
                    elseif II1==size(W4d,1)
                        delta_4_gs_ms(:,2,COUNT)=[(Momentstable(1:9,II1-1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))/(Momentstable(8,II1-1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(8,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))];
                    elseif II1==1
                        delta_4_gs_ms(:,2,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1+1,deltaset(II2),rhoset(II3),vssset(II4),II5))/(Momentstable(8,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(8,II1+1,deltaset(II2),rhoset(II3),vssset(II4),II5))];
                    end
                    if size(W4d,5)>1
                        if II5<size(W4d,5) && II5>1
                            delta_4_gs_ms(:,3,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5-1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5+1))/(Momentstable(8,II1,deltaset(II2),rhoset(II3),vssset(II4),II5-1)-Momentstable(8,II1,deltaset(II2),rhoset(II3),vssset(II4),II5+1))];
                        elseif II5==size(W4d,5) && size(W4d,5) >1
                            delta_4_gs_ms(:,3,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5-1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))/(Momentstable(8,II1,deltaset(II2),rhoset(II3),vssset(II4),II5-1)-Momentstable(8,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))];
                        elseif II5==1
                            delta_4_gs_ms(:,3,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5+1))/(Momentstable(8,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(8,II1,deltaset(II2),rhoset(II3),vssset(II4),II5+1))];
                        end
                        Vdelta_4_gs_ms(:,COUNT)=diag(inv((delta_4_gs_ms(4:7,:,COUNT)'*inv((VCV2data(1:4,1:4)) ))*delta_4_gs_ms(4:7,:,COUNT)));
                    else
                        Vdelta_4_gs_ms(:,COUNT)=diag(inv((delta_4_gs_ms(4:7,1:2,COUNT)'*inv((VCV2data(1:4,1:4)) ))*delta_4_gs_ms(4:7,1:2,COUNT)));
                    end
                    
                    % ===============
                    % two moments
                    % ===============
                    WW=W2d(:,deltaset,rhoset,vssset,:);
                    criterion=reshape(WW,prod(size(WW)),1);
                    [CC,II]=min(criterion);
                    [II1,II2,II3,II4,II5]=ind2sub(size(WW),II);
                    solution_2_gs_ms(:,COUNT)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5); CC];
                    if II2<length(deltaset) && II2>1
                        delta_2_gs_ms(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))];
                    elseif II2==length(deltaset)
                        delta_2_gs_ms(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))];
                    elseif II2==1
                        delta_2_gs_ms(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))];
                    end
                    
                    Vdelta_2_gs_ms(:,COUNT)=inv(delta_2_gs_ms(4:5,COUNT)'*inv((VCV2data(1:2,1:2)))*delta_2_gs_ms(4:5,COUNT));
                    
                    % ===============
                    % find minimum over all groupsizes, no measurement error
                    % ===============
                    % four moments
                    WW=W4d(1,deltaset,rhoset,vssset,1);
                    criterion=reshape(WW,prod(size(WW)),1);
                    [CC,II]=min(criterion);
                    [II1,II2,II3,II4,II5]=ind2sub(size(WW),II);
                    solution_4_gs(:,COUNT)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5); CC];
                    if II2<length(deltaset) && II2>1
                        delta_4_gs(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))];
                    elseif II2==length(deltaset)
                        delta_4_gs(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))];
                    elseif II2==1
                        delta_4_gs(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))];
                    end
                    Vdelta_4_gs(:,COUNT)=diag(inv(delta_4_gs(4:7,COUNT)'*inv((VCV2data(1:4,1:4)))*delta_4_gs(4:7,COUNT)));
                    
                    % two moments
                    WW=W2d(1,deltaset,rhoset,vssset,1);
                    criterion=reshape(WW,prod(size(WW)),1);
                    [CC,II]=min(criterion);
                    [II1,II2,II3,II4,II5]=ind2sub(size(WW),II);
                    solution_2_gs(:,COUNT)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5); CC];
                    if II2<length(deltaset) && II2>1
                        delta_2_gs(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))];
                    elseif II2==length(deltaset)
                        delta_2_gs(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vssset(II4),II5))];
                    elseif II2==1
                        delta_2_gs(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vssset(II4),II5)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vssset(II4),II5))];
                    end
                    Vdelta_2_gs(:,COUNT)=inv(delta_2_gs(4:5,COUNT)'*inv((VCV2data(1:2,1:2)))*delta_2_gs(4:5,COUNT));
                    
                    
                    % same for maximum group size
                    if coaldev==COALDEV(1)
                        delta_4_ms=nan(9,3,2);delta_4=nan(9,1,2);delta_2_ms=nan(9,3,2);delta_2=nan(9,1,2);
                    end
                    
                    % four moments with measurement error
                    WW=W4dm(:,deltaset,rhoset,:);
                    criterion=reshape(WW,prod(size(WW)),1);
                    [CC,II]=min(criterion);
                    [II1,II2,II3,II4]=ind2sub(size(WW),II);
                    vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)).*(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)>0)));
                    solution_4_ms(:,COUNT)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4) ; CC];
                    if II2<length(deltaset) && II2>1 && isnan(Momentstable(1,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))==0
                        delta_4_ms(:,1,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                    elseif II2==length(deltaset)
                        delta_4_ms(:,1,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                    elseif II2==1 || isnan(Momentstable(1,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))==1
                        delta_4_ms(:,1,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                    end
                    
                    if II1<size(W4dm,1) && II1>1
                        delta_4_ms(:,2,COUNT)=[(Momentstable(1:9,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(8,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(8,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))];
                    elseif II1==size(W4dm,1)
                        delta_4_ms(:,2,COUNT)=[(Momentstable(1:9,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(8,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(8,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                    elseif II1==1
                        delta_4_ms(:,2,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(8,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(8,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))];
                    end
                    if size(W4dm,4)>1
                        if II4<size(W4dm,4) && II4>1
                            delta_4_ms(:,3,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))];
                        elseif II4==size(W4dm,4)
                            delta_4_ms(:,3,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                        elseif II4==1
                            delta_4_ms(:,3,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))];
                        end
                    end
                    
                    if sum(isnan(delta_4_ms(:,3,COUNT)))==0
                        
                        Vdelta_4_ms(:,COUNT)=diag(inv((delta_4_ms(4:7,:,COUNT)'*inv((VCV2data(1:4,1:4)) ))*delta_4_ms(4:7,:,COUNT)));
                    else
                        Vdelta_4_ms(:,COUNT)=[diag(inv((delta_4_ms(4:7,1:2,COUNT)'*inv((VCV2data(1:4,1:4)) ))*delta_4_ms(4:7,1:2,COUNT)));nan];
                    end
                    
                    % indicator for estimated meas error
                    III4(village,coaldev+1)=II4;
                    III1(village,coaldev+1)=II1;
                    III2(village,coaldev+1)=II2;

                    % measurement error only in consumption
                    solution_4_ms_plus(:,COUNT+Countadd)=solution_4_ms(:,COUNT);
                    delta_4_ms_plus(:,3,COUNT+Countadd)=delta_4_ms(:,3,COUNT);
                    Vdelta_4_ms_plus(:,COUNT+Countadd)=Vdelta_4_ms(:,COUNT);
                    if coaldev==1
                        Countadd=Countadd+1;
                        WW=W4dm(:,deltaset,rhoset,1);
                        criterion=reshape(WW,prod(size(WW)),1);
                        [CC,II]=min(criterion);
                        [II1,II2,II3]=ind2sub(size(WW),II); II4=1;
                        vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)).*(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)>0)));
                        solution_4_ms_plus(:,COUNT+Countadd)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4) ; CC];
                        if II2<length(deltaset) && II2>1 && isnan(Momentstable(1,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))==0
                            delta_4_ms_plus(:,1,COUNT+Countadd)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                        elseif II2==length(deltaset)
                            delta_4_ms_plus(:,1,COUNT+Countadd)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                        elseif II2==1 || isnan(Momentstable(1,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))==1
                            delta_4_ms_plus(:,1,COUNT+Countadd)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                        end
                        
                        if II1<size(W4dm,1) && II1>1
                            delta_4_ms_plus(:,2,COUNT+Countadd)=[(Momentstable(1:9,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(8,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(8,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))];
                        elseif II1==size(W4dm,1)
                            delta_4_ms_plus(:,2,COUNT+Countadd)=[(Momentstable(1:9,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(8,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(8,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                        elseif II1==1
                            delta_4_ms_plus(:,2,COUNT+Countadd)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(8,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(8,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))];
                        end
                        if size(W4dm,4)>1
                            if II4<size(W4dm,4) && II4>1
                                delta_4_ms_plus(:,3,COUNT+Countadd)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))];
                            elseif II4==size(W4dm,4)
                                delta_4_ms_plus(:,3,COUNT+Countadd)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                            elseif II4==1
                                delta_4_ms_plus(:,3,COUNT+Countadd)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))];
                            end
                        end
                        
                        if sum(isnan(delta_4_ms_plus(:,3,COUNT+Countadd)))==0
                            Vdelta_4_ms_plus(:,COUNT+Countadd)=diag(inv((delta_4_ms_plus(4:7,:,COUNT+Countadd)'*inv((VCV2data(1:4,1:4)) ))*delta_4_ms_plus(4:7,:,COUNT+Countadd)));
                        else
                            Vdelta_4_ms_plus(:,COUNT+Countadd)=[diag(inv((delta_4_ms_plus(4:7,1:2,COUNT+Countadd)'*inv((VCV2data(1:4,1:4)) ))*delta_4_ms_plus(4:7,1:2,COUNT+Countadd)));nan];
                        end
                    end
                    
                    % two moments with measurement error
                    WW=W2dm(:,deltaset,rhoset,:);
                    criterion=reshape(WW,prod(size(WW)),1);
                    [CC,II]=min(criterion);
                    [II1,II2,II3,II4]=ind2sub(size(WW),II);
                    vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4).*(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)>0))));
                    solution_2_ms(:,COUNT)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4) ; CC];
                    if II2<length(deltaset) && II2>1
                        delta_2_ms(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                    elseif II2==length(deltaset)
                        delta_2_ms(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                    elseif II2==1
                        delta_2_ms(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                    end
                    Vdelta_2_ms(:,COUNT)=inv(delta_2_ms(4:5,COUNT)'*inv((VCV2data(1:2,1:2)))*delta_2_ms(4:5,COUNT));
                end
                
                if est_rho1==1 % for estimated persistence
                    % four moments with measurement error
                    WW=W4dm(1,deltaset,rhoset,:);
                    criterion=reshape(WW,prod(size(WW)),1);
                    [CC,II]=min(criterion);
                    [II1,II2,II3,II4]=ind2sub(size(WW),II);
                    vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)).*(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)>0)));
                    solution_4_rho1(:,COUNT)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4) ; CC];
                    
                    % NB: here the change in delta often implies changes in sustainable
                    % group size
                    if II2<length(deltaset) && II2>1
                        if isnan(Momentstable(1,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))==0
                            delta_4_rho1(:,1,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                        else
                            delta_4_rho1(:,1,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                        end
                    elseif II2==length(deltaset)
                        if isnan(Momentstable(1,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))==0
                            delta_4_rho1(:,1,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                        else
                            delta_4_rho1(:,1,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max-1,II4)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max-1,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max-1,II4)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max-1,II4))];
                        end
                    elseif II2==1
                        delta_4_rho1(:,1,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                    end
                    
                    if II4<size(W4dm,4) && II4>1
                        if isnan(Momentstable(4,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))==0
                            delta_4_rho1(:,2,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))];
                        else
                            delta_4_rho1(:,2,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                        end
                    elseif II4==size(W4dm,4)
                        delta_4_rho1(:,2,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                    elseif II4==1 || isnan(Momentstable(4,II1,deltaset(II2),rhoset(II3),vss_max,II4-1))==1
                        delta_4_rho1(:,2,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))];
                    end
                    Vdelta_4_rho1(:,COUNT)=diag(inv((delta_4_rho1(4:7,:,COUNT)'*inv((VCV2data(1:4,1:4)) ))*delta_4_rho1(4:7,:,COUNT)))
                end
                
                % four moments no measurement error
                if coaldev==0 && buildup==1 && samedelta==1
                    II2=deltaset(	squeeze(Momentstable(2,1,:,1,max(vss_set),1))==solution_4(2,COUNT-1));
                    II3=1; CC=nan(length(vss_set),1)'; vss_max=vss_set;
                    Vdelta_4(:,COUNT:COUNT+length(buildupvec)-1)=nan(length(Vdelta_4(:,COUNT-1)),length(vss_set));
                    solution_4(:,COUNT:COUNT+length(buildupvec)-1)=[squeeze(Momentstable(1:9,1,deltaset(II2),rhoset(II3),vss_max,1) ); CC];
                else
                    WW=W4dm(1,deltaset,rhoset,1);
                    criterion=reshape(WW,prod(size(WW)),1);
                    [CC,II]=min(criterion);
                    [II1,II2,II3,II4]=ind2sub(size(WW),II);
                    vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,1))));
                    if II2>1
                        vss_max1=max(find(isfinite(Momentstable(1,II1,deltaset(II2-1),rhoset(II3),:,1))));
                    end
                    if II2<length(deltaset)
                        vss_max2=max(find(isfinite(Momentstable(1,II1,deltaset(II2+1),rhoset(II3),:,1))));
                    end
                    
                    solution_4(:,COUNT)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,1) ; CC];
                    if II2<length(deltaset) && II2>1
                        delta_4(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))];
                    elseif II2==length(deltaset)
                        delta_4(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,1))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,1))];
                    elseif II2==1
                        delta_4(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,1)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,1)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))];
                    end
                    Vdelta_4(:,COUNT)=inv(delta_4(4:7,COUNT)'*inv((VCV2data(1:4,1:4)))*delta_4(4:7,COUNT));
                end
                
                % two moments no measurement error
                WW=W2dm(1,deltaset,rhoset,1);
                criterion=reshape(WW,prod(size(WW)),1);
                [CC,II]=min(criterion);
                [II1,II2,II3]=ind2sub(size(WW),II);
                vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,1))));
                if II2>1
                    vss_max1=max(find(isfinite(Momentstable(1,II1,deltaset(II2-1),rhoset(II3),:,1))));
                end
                
                if II2<length(deltaset)
                    vss_max2=max(find(isfinite(Momentstable(1,II1,deltaset(II2+1),rhoset(II3),:,1))));
                end
                solution_2(:,COUNT)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,1) ; CC];
                if II2<length(deltaset) && II2>1
                    delta_2(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))];
                elseif II2==length(deltaset)
                    delta_2(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,1))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,1))];
                elseif II2==1
                    delta_2(:,COUNT)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,1)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,1)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))];
                end
                Vdelta_2(:,COUNT)=inv(delta_2(4:5,COUNT)'*inv(diag(diag(VCV2data(1:2,1:2))))*delta_2(4:5,COUNT));
                if coaldev==1 && village==1
                    %stop
                end
                %find minimum for village/largest stable group, no
                %measurement error, raw moments
                % find for group size==3, no measurement error ,
                %%%four moments
                
                WW=W4d(1,:,:,2,1);
                criterion=reshape(WW,prod(size(WW)),1);
                [CC,II]=min(criterion);
                [II1,II2,II3,II4,II5]=ind2sub(size(WW),II);
                solution_4_gs3(:,COUNT)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),3,II5); CC];
                
            elseif coaldev==2 % this is for the savings economy
                delta_4_ms_SI=nan(10,4,1);
                WW=W4dm(:,deltaset,rhoset,:);
                WW=WW.*(WW./WW); % eliminate zeros
                criterion=reshape(WW,prod(size(WW)),1);
                [CC,II]=min(criterion);
                [II1,II2,II3,II4]=ind2sub(size(WW),II);
                vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)).*(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)>0)));
                solution_4_ms_SI(:,1)=[Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4) ; CC];
                if II2<length(deltaset) && II2>1
                    delta_4_ms_SI(:,1)=[(Momentstable(:,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                elseif II2==length(deltaset)
                    delta_4_ms_SI(:,1)=[(Momentstable(:,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                elseif II2==1
                    delta_4_ms_SI(:,1)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                end
                delta_4_ms_SI(1,1)=nan;
                if II1<size(W4dm,1) && II1>1
                    delta_4_ms_SI(:,2)=[(Momentstable(:,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(8,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(8,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))];
                elseif II1==size(W4dm,1)
                    delta_4_ms_SI(:,2)=[(Momentstable(:,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(8,II1-1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(8,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                elseif II1==1
                    delta_4_ms_SI(:,2)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(8,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(8,II1+1,deltaset(II2),rhoset(II3),vss_max,II4))];
                end
                
                if size(W4dm,4)>1
                    if II4<size(W4dm,4) && II4>1
                        delta_4_ms_SI(:,3)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))];
                    elseif II4==size(W4dm,4)
                        delta_4_ms_SI(:,3)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4-1)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                    elseif II4==1
                        delta_4_ms_SI(:,3)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))/(Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(9,II1,deltaset(II2),rhoset(II3),vss_max,II4+1))];
                    end
                else
                    delta_4_ms_SI(:,3)=nan(length(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4)),1);
                end
                
                % indicator for estimated meas error
                III4(village,coaldev+1)=II4;
                III1(village,coaldev+1)=II1;
                III2(village,coaldev+1)=II2;
                if length(rhoset)==1
                    
                elseif II3<length(rhoset) && II3>1
                    delta_4_ms_SI(:,4)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))/(Momentstable(10,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(10,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))];
                elseif II3==length(rhoset) && II3>0
                    delta_4_ms_SI(:,4)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(10,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(10,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                elseif II3==1
                    delta_4_ms_SI(:,4)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))/(Momentstable(10,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(10,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))];
                end
                
                Vdelta_4_ms_SI(:)=diag(inv((delta_4_ms_SI(4:7,~isnan(delta_4_ms_SI(5,:)))'*inv((VCV2data(1:4,1:4)) ))*delta_4_ms_SI(4:7,~isnan(delta_4_ms_SI(5,:)))));
                clear WW
                WW=W4dm(1,deltaset,rhoset,1);
                WW=WW.*(WW./WW); % eliminate zeros
                criterion=reshape(WW,prod(size(WW)),1);
                [CC,II]=min(criterion);
                [II1,II2,II3,II4]=ind2sub(size(WW),II);
                vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)).*(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)>0)));
                solution_4_SI(:,1)=[Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4) ; CC];
                
                if II2<length(deltaset) && II2>1
                    delta_4_SI(:,1)=[(Momentstable(:,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                elseif II2==length(deltaset)
                    delta_4_SI(:,1)=[(Momentstable(:,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                elseif II2==1
                    delta_4_SI(:,1)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                end
                delta_4_SI(1,1)=nan;
                
                if length(rhoset)==1
                elseif II3<length(rhoset) && II3>1
                    delta_4_SI(:,2)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))/(Momentstable(10,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(10,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))];
                elseif II3==length(rhoset) && II3>0
                    delta_4_SI(:,2)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(10,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(10,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                elseif II3==1
                    delta_4_SI(:,2)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))/(Momentstable(10,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(10,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))];
                end
                Vdelta_4_SI(:)=diag(inv((delta_4_SI(4:7,:)'*inv((VCV2data(1:4,1:4)) ))*delta_4_SI(4:7,:)));
                WW=W2dm(1,deltaset,rhoset,1);
                WW=WW.*(WW./WW); % eliminate zeros
                criterion=reshape(WW,prod(size(WW)),1);
                [CC,II]=min(criterion);
                [II1,II2,II3,II4]=ind2sub(size(WW),II);
                vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)).*(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)>0)));
                solution_2_SI(:,1)=[Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4) ; CC];
                
                if II2<length(deltaset) && II2>1
                    delta_2_SI(:,1)=[(Momentstable(:,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                elseif II2==length(deltaset)
                    delta_2_SI(:,1)=[(Momentstable(:,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                elseif II2==1
                    delta_2_SI(:,1)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                end
                delta_2_SI(1,1)=nan;
                
                if length(rhoset)==1
                elseif II3<length(rhoset) && II3>1
                    delta_2_SI(:,2)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))/(Momentstable(10,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(10,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))];
                elseif II3==length(rhoset) && II3>0
                    delta_2_SI(:,2)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(10,II1,deltaset(II2),rhoset(II3-1),vss_max,II4)-Momentstable(10,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                elseif II3==1
                    delta_2_SI(:,2)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))/(Momentstable(10,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(10,II1,deltaset(II2),rhoset(II3+1),vss_max,II4))];
                end
                
                Vdelta_2_SI(:)=diag(inv((delta_2_SI(4:5,:)'*inv((VCV2data(1:2,1:2)) ))*delta_2_SI(4:5,:)));
                
            elseif coaldev==3
                clear vss_max_j
                WW=W4dm(1,deltaset,rhoset,1);
                criterion=reshape(WW,prod(size(WW)),1);
                [CC,II]=min(criterion);
                [II1,II2,II3,II4]=ind2sub(size(WW),II);
                vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)).*(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)>0)));
                for jjj=1:size(Momentstable,3)
                    if isempty(max(find(isfinite(Momentstable(1,II1,jjj,rhoset(II3),:,II4)).*(Momentstable(1,II1,jjj,rhoset(II3),:,II4)>0))))==0
                        vss_max_j(jjj)=max(find(isfinite(Momentstable(1,II1,jjj,rhoset(II3),:,II4)).*(Momentstable(1,II1,jjj,rhoset(II3),:,II4)>0)));
                    else
                        vss_max_j(jjj)=0;
                    end
                end
                solution_4_FI(:,1)=[Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4) ; CC];
                
                delta1temp=min(squeeze(Momentstable(2,II1,find(vss_max_j==vss_max),rhoset(II3),vss_max,II4)));
                delta2temp=max(squeeze(Momentstable(2,II1,find(vss_max_j==vss_max),rhoset(II3),vss_max,II4)));
                
            elseif coaldev==4
                
                WW=W4dm(1,deltaset,rhoset,1);
                criterion=reshape(WW,prod(size(WW)),1);
                [CC,II]=min(criterion);
                [II1,II2,II3,II4]=ind2sub(size(WW),II);
                vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)).*(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,II4)>0)));
                solution_4_het(:,1)=[Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4) ; CC];
                
                if length(deltaset)>1 && II2<length(deltaset) && II2>1
                    if isnan(Momentstable(1,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))==0 && isnan(Momentstable(1,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))==0
                        delta_4_het(:,1)=[(Momentstable(:,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))];
                    elseif isnan(Momentstable(1,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))==0
                        delta_4_het(:,1)=[(Momentstable(:,II1,deltaset(II2+1),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                    else
                        delta_4_het(:,1)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))];
                    end
                elseif length(deltaset)>1 &&  II2==length(deltaset) && isnan(Momentstable(1,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))==0
                    delta_4_het(:,1)=[(Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max,II4))];
                elseif length(deltaset)>1 && II2==1 &&  isnan(Momentstable(1,II1,deltaset(II2+1),rhoset(II3),vss_max,II4))==0
                    delta_4_het(:,1)=[(Momentstable(:,II1,deltaset(II2+1),rhoset(II3),vss_max,II4)-Momentstable(:,II1,deltaset(II2),rhoset(II3),vss_max,II4))/(Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max,II4)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,II4))];
                else
                    delta_4_het(:,1)=nan;
                end
                if      isnan(delta_4_het(1,1))==0
                    Vdelta_4_het(:)=diag(inv(delta_4_het(4:7,find(~isnan(delta_4_het(5,:))))'*inv((VCV2data(1:4,1:4)) )*delta_4_het(4:7,find(~isnan(delta_4_het(5,:))))));
                else
                    Vdelta_4_het(:)=[nan];
                end
                
                % two moments no measurement error
                WW=W2dm(1,deltaset,rhoset,1);
                criterion=reshape(WW,prod(size(WW)),1);
                [CC,II]=min(criterion);
                [II1,II2,II3]=ind2sub(size(WW),II);
                vss_max=max(find(isfinite(Momentstable(1,II1,deltaset(II2),rhoset(II3),:,1))));
                if II2>1
                    vss_max1=max(find(isfinite(Momentstable(1,II1,deltaset(II2-1),rhoset(II3),:,1))));
                end
                
                if II2<length(deltaset)
                    vss_max2=max(find(isfinite(Momentstable(1,II1,deltaset(II2+1),rhoset(II3),:,1))));
                end
                solution_2_het(:,1)=[Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,1) ; CC];
                if II2<length(deltaset) && II2>1
                    delta_2_het(:,1)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))];
                elseif II2==length(deltaset)
                    delta_2_het(:,1)=[(Momentstable(1:9,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,1))/(Momentstable(2,II1,deltaset(II2-1),rhoset(II3),vss_max1,1)-Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,1))];
                elseif II2==1
                    delta_2_het(:,1)=[(Momentstable(1:9,II1,deltaset(II2),rhoset(II3),vss_max,1)-Momentstable(1:9,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))/(Momentstable(2,II1,deltaset(II2),rhoset(II3),vss_max,1)-Momentstable(2,II1,deltaset(II2+1),rhoset(II3),vss_max2,1))];
                end
                Vdelta_2_het(:,1)=inv(delta_2_het(4:5,1)'*inv(diag(diag(VCV2data(1:2,1:2))))*delta_2_het(4:5,1));
            end
        end
    end
    Momentsdatatable=[nan;nan;nan;Momvec2data'];
    
    if rhogrid==1
        disp('only do identification for 1 village - so stop here and run identif_pic')
        stop
    else
        if sum(COALDEV-1==0)==1 && sum(COALDEV==0)==0
            Moments_4_table(:,(village-1)*2+1:village*2)=[[nan;Momentsdatatable;nan],[solution_4([1:2],:);Vdelta_4.^0.5;solution_4([3:7, 10],:)]];
        elseif sum(COALDEV-1==0)==1 && sum(COALDEV==0)==1
            if sum(rhovec==rho_calib)==length(rhovec)
                if relvarincscale==1
                    Momentscalibtable(:,(village-1)*(size(Momentscalib,2)+1)+1:village*(size(Momentscalib,2)+1))=[Momentsdatatable,[Momentscalib([1,2,3,8,15],:);Momentscalib(9,:)/Vardylog_cond(village);Momentscalib(16,:)-Momentscalib(17,:)]];
                else
                    Momentscalibtable(:,(village-1)*(size(Momentscalib,2)+1)+1:village*(Momentscalib(COALDEV,2)+1))=[Momentsdatatable,[Moments([1,2,3,8,15],:);Momentscalib(9,:)./Momentscalib(8,:)./Vardylog_cond(village);Momentscalib(16,:)-Momentscalib(17,:)]];
                end
            end
            Moments_2_table(:,(village-1)*3+1:village*3)=[[nan;Momentsdatatable;nan],[solution_2([1:2],:);Vdelta_2.^0.5;solution_2([3:7, 10],:)]];
            % Moments_4_table(:,(village-1)*3+1:village*3)=[[Momentsdatatable;nan],solution_4([1:7, 10],:)];
            if buildup==0
                Moments_4_table(:,(village-1)*3+1:village*3)=[[nan;Momentsdatatable;nan],[solution_4([1:2],:);Vdelta_4.^0.5;solution_4([3:7, 10],:)]];
            elseif buildup==1
                Moments_4_table(:,1:length(buildupvec)+2)=[[nan;Momentsdatatable;nan],[solution_4([1:2],:);Vdelta_4.^0.5;solution_4([3:7, 10],:)]];
            end
            %     Moments_2_res_table(:,(village-1)*3+1:village*3)=[[Momentsdatatable;nan],solution_2_res([1:7, 10],:)];
            %     Moments_4_res_table(:,(village-1)*3+1:village*3)=[[Momentsdatatable;nan],solution_4_res([1:7, 10],:)];
            
            if est_rho1==1
                Moments_4_rho1_table(:,(village-1)*3+1:village*3)=[[nan;Momentsdatatable;nan;nan;nan],[solution_4_rho1([1:2],:);Vdelta_4_rho1(1,:).^0.5;solution_4_rho1([3:7],:);solution_4_rho1([9],:);Vdelta_4_rho1(2,:).^0.5;solution_4_rho1([10],:)]];
            end
            if Nincproc>1
                Moments_4_ms_table(:,(village-1)*3+1:village*3)=[[nan;Momentsdatatable;nan;nan;nan;nan;nan],[solution_4_ms([1:2],:);Vdelta_4_ms(1,:).^0.5;solution_4_ms([3:8],:);Vdelta_4_ms(2,:).^0.5;solution_4_ms([9],:);Vdelta_4_ms(3,:).^0.5;solution_4_ms([10],:)]];
                Moments_4_gs_table(:,(village-1)*3+1:village*3)=[[Momentsdatatable;nan],solution_4_gs([1:7, 10],:)];
                Moments_4_gs3_table(:,(village-1)*3+1:village*3)=[[Momentsdatatable;nan],[solution_4_gs3([1:7, 10],:)]]; % solution_4([1:7, 10],2)
            end
        end
        
        if sum(COALDEV-1==0)==0 && sum(COALDEV==0)==1
            if est_rho1==1
                Moments_4_rho1_table(:,(village-1)*2+1:village*2)=[[nan;Momentsdatatable;nan;nan;nan],[solution_4_rho1([1:2],:);Vdelta_4_rho1(1,:).^0.5;solution_4_rho1([3:7],:);solution_4_rho1([9],:);Vdelta_4_rho1(2,:).^0.5;solution_4_rho1([10],:)]];
            end
        end
        
        if sum(COALDEV-2==0)==1
            if ext==1
                if isnan(phicalib)==1
                    Moments_4_table_plus(:,(village-1)*4+1:village*4)=[[nan;nan;nan;Momentsdatatable;nan],[solution_4([1:2],:);Vdelta_4.^0.5;solution_4([3],:);[nan,nan;nan,nan];solution_4([4:7, 10],:)],[nan;solution_4_SI([2],:);Vdelta_4_SI(1).^0.5;solution_4_SI([3,10],:);Vdelta_4_SI(2).^0.5;solution_4_SI([4:7, 11],:)]];
                    Moments_2_table_plus(:,(village-1)*4+1:village*4)=[[nan;nan;nan;Momentsdatatable;nan],[solution_2([1:2],:);Vdelta_2.^0.5;solution_4([3],:);[nan,nan;nan,nan];solution_2([4:7, 10],:)],[nan;solution_2_SI([2],:);Vdelta_2_SI(1).^0.5;solution_2_SI([3,10],:);Vdelta_2_SI(2).^0.5;solution_2_SI([4:7, 11],:)]];
                    Moments_4_ms_table_plus(:,(village-1)*4+1:village*4)=[[nan;nan;nan;Momentsdatatable;nan;nan;nan;nan;nan],[solution_4_ms([1:2],:);Vdelta_4_ms(1,:).^0.5;solution_4_ms([3],:);[nan,nan;nan,nan];solution_4_ms([4:8],:);Vdelta_4_ms(2,:).^0.5;solution_4_ms(9,:);Vdelta_4_ms(3,:).^0.5;solution_4_ms([10],:)],[nan;solution_4_ms_SI([2],:);Vdelta_4_ms_SI(1).^0.5;solution_4_ms_SI([3,10],:);Vdelta_4_ms_SI(2).^0.5;solution_4_ms_SI([4:8],:);Vdelta_4_ms_SI(3).^0.5;solution_4_ms_SI(9,:);Vdelta_4_ms_SI(4).^0.5;solution_4_ms_SI([11],:)]];
                else
                    Moments_4_table_plus(:,(village-1)*4+1:village*4)=[[nan;Momentsdatatable;nan],[solution_4([1:2],:);Vdelta_4.^0.5;solution_4([3],:);solution_4([4:7, 10],:)],[nan;solution_4_SI([2],:);Vdelta_4_SI(1).^0.5;solution_2_SI([3],:);solution_4_SI([4:7, 11],:)]];
                    Moments_2_table_plus(:,(village-1)*4+1:village*4)=[[nan;Momentsdatatable;nan],[solution_2([1:2],:);Vdelta_2.^0.5;solution_4([3],:);solution_2([4:7, 10],:)],[nan;solution_2_SI([2],:);Vdelta_2_SI(1).^0.5;solution_2_SI([3],:);solution_2_SI([4:7, 11],:)]];
                    
                    if Nincproc>1
                        Moments_4_ms_table_plus(:,(village-1)*4+1:village*4)=[[nan;nan;Momentsdatatable;nan;nan;nan;nan;nan],[solution_4_ms([1:2],:);Vdelta_4_ms(1,:).^0.5;solution_4_ms([3],:);[nan,nan];solution_4_ms([4:8],:);Vdelta_4_ms(2,:).^0.5;solution_4_ms(9,:);Vdelta_4_ms(3,:).^0.5;solution_4_ms([10],:)],[nan;solution_4_ms_SI([2],:);Vdelta_4_ms_SI(1).^0.5;solution_4_ms_SI([3,10],:);solution_4_ms_SI([4:8],:);Vdelta_4_ms_SI(2).^0.5;solution_4_ms_SI(9,:);Vdelta_4_ms_SI(3).^0.5;solution_4_ms_SI([11],:)]];
                        Moments_4_ms_table_plus_plus(:,(village-1)*5+1:village*5)=[[nan;nan;Momentsdatatable;nan;nan;nan;nan;nan],[solution_4_ms_plus([1:2],:);Vdelta_4_ms_plus(1,:).^0.5;solution_4_ms_plus([3],:);[nan,nan,nan];solution_4_ms_plus([4:8],:);Vdelta_4_ms_plus(2,:).^0.5;solution_4_ms_plus(9,:);Vdelta_4_ms_plus(3,:).^0.5;solution_4_ms_plus([10],:)],[nan;solution_4_ms_SI([2],:);Vdelta_4_ms_SI(1).^0.5;solution_4_ms_SI([3,10],:);solution_4_ms_SI([4:8],:);Vdelta_4_ms_SI(2).^0.5;solution_4_ms_SI(9,:);Vdelta_4_ms_SI(3).^0.5;solution_4_ms_SI([11],:)]];
                    end
                end               
            end
            if isnan(phicalib)==1
                Moments_4_ms_table_SI(:,(village-1)*2+1:village*2)=[[nan;nan;nan;Momentsdatatable;nan;nan;nan;nan;nan],[solution_4_ms_SI([1:2],:);Vdelta_4_ms_SI(1).^0.5;solution_4_ms_SI(3,:);solution_4_ms_SI(10,:);Vdelta_4_ms_SI(4).^0.5;solution_4_ms_SI([4:8],:);Vdelta_4_ms_SI(2).^0.5;solution_4_ms_SI([9],:);Vdelta_4_ms_SI(3).^0.5;solution_4_ms_SI([11],:)]];;
                Moments_4_table_SI(:,(village-1)*2+1:village*2)=[[nan;nan;nan;Momentsdatatable;nan],[solution_4_SI([1:2],:);Vdelta_4_SI(1).^0.5;solution_4_ms_SI(3,:);;solution_4_SI(10,:);Vdelta_4_ms_SI(2).^0.5;solution_4_SI([4:7],:);solution_4_ms_SI([11],:)]];
                Moments_2_table_SI(:,(village-1)*2+1:village*2)=[[nan;nan;nan;Momentsdatatable;nan],[solution_2_SI([1:2],:);Vdelta_2_SI(1).^0.5;solution_4_ms_SI(3,:);;solution_2_SI(10,:);Vdelta_4_ms_SI(2).^0.5;solution_2_SI([4:7],:);solution_4_ms_SI([11],:)]];
            end
            
        end
        if sum(COALDEV-3==0)==1
            %  Momentscalibtable_FI(:,(village-1)*2+1:village*2)=[Momentsdatatable,[Momentscalib([1,2,3,8,15],COUNT);Momentscalib(9,COUNT)/Vardylog_cond(village);Momentscalib(16,COUNT)-Momentscalib(17,COUNT)]];
            
            Moments_4_table_FI(:,(village-1)*2+1:village*2)=[[nan;Momentsdatatable;nan],[solution_4_FI([1],:);delta1temp;delta2temp;solution_4_FI(3,:);solution_4_FI([4:7],:);solution_4_FI([11],:)]];;
        end
        
        if sum(COALDEV==0)==1 && sum(COALDEV-4==0)==1
            Moments_4_table_het(:,(village-1)*3+1:village*3)=[[nan;Momentsdatatable;nan],[solution_4([1:2],1);Vdelta_4(1).^0.5;solution_4([3:7, 10],1)],[solution_4_het([1:2],1);Vdelta_4_het(1).^0.5;solution_4_het([3:7, 11],1)]]
            Moments_2_table_het(:,(village-1)*3+1:village*3)=[[nan;Momentsdatatable;nan],[solution_2([1:2],1);Vdelta_2(1).^0.5;solution_2([3:7, 10],1)],[solution_2_het([1:2],1);Vdelta_2_het(1).^0.5;solution_2_het([3:7, 10],1)]]
            Moments_4_table_het=real(Moments_4_table_het);
        end
    end
end

if write==1 && rhogrid==0
    if Tessa1Tobi2==1
        path1= [ currentfolder '\Paper\Latex\'];
    else
        path1= [ currentfolder '\Paper\Latex\'];
    end
    % pre-chosen parameters
    if sum(COALDEV-1==0)==1 && sum(COALDEV==0)==1
        if buildup==1
            if COALDEV(1)==0
                columnLabels = {'Data','ID', 'CD','Data','ID', 'CD','Data','ID', 'CD'};
            else
                columnLabels = {'Data', 'CD','ID','ID','ID'};
            end
            
            % parameters chosen to target 2 moments
            if pres==1
                rowLabels = {'n','$\delta$','s.e.','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\Delta V_{dc}$', '$\Delta \beta_{dc}$', 'Fit'};
            else
                rowLabels = {'n','$\delta$','s.e.','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy\leq0}$', 'Model criterion'};
            end
            matrix2latex_newnew(Moments_2_table,[path1 'same_vss_moments2_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{3}{c}{\textbf{Aurepalle}}& \multicolumn{3}{c}{\textbf{Kanzara}}& \multicolumn{3}{c}{\textbf{Shirapur}}');
            % parameters chosen to target 4 moments
            matrix2latex_newnew(Moments_4_table,[path1 'same_vss_moments4_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{3}{c}{\textbf{Aurepalle}}');
        else
            if sum(rhovec==rho_calib)==length(rhovec)
                rowLabels = {'n','$\delta$','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy \leq 0}$'};
                if COALDEV(1)==0
                    columnLabels = {'Data','ID', 'CD','SI','Data','ID', 'CD','SI','Data','ID', 'CD','SI'};
                else
                    columnLabels = {'Data', 'CD','ID','SI','Data','CD','ID','SI','Data','CD','ID','SI'};
                end
                path=[path1 'momentscalib_' int2str(outlier) int2str(fixedeffect) Date '.txt'];
                matrix2latex_newnew(Momentscalibtable,path,'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{3}{c}{\textbf{Aurepalle}}& \multicolumn{3}{c}{\textbf{Kanzara}}& \multicolumn{3}{c}{\textbf{Shirapur}}');
                
            end
            if COALDEV(1)==0
                columnLabels = {'Data','ID', 'CD','Data','ID', 'CD','Data','ID', 'CD'};
            else
                columnLabels = {'Data', 'CD','ID','Data','CD','ID','Data','CD','ID'};
            end
            
            % parameters chosen to target 2 moments
            if pres==1
                rowLabels = {'n','$\delta$','s.e.','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\Delta V_{dc}$', '$\Delta \beta_{dc}$', 'Fit'};
            else
                rowLabels = {'n','$\delta$','s.e.','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy\leq0}$', 'Model criterion'};
            end
            matrix2latex_newnew(Moments_2_table,[path1 'moments2_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{3}{c}{\textbf{Aurepalle}}& \multicolumn{3}{c}{\textbf{Kanzara}}& \multicolumn{3}{c}{\textbf{Shirapur}}');
            % parameters chosen to target 4 moments
            if VILLAGE==[1,2,3]
                matrix2latex_newnew(Moments_4_table,[path1 'moments4_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{3}{c}{\textbf{Aurepalle}}& \multicolumn{3}{c}{\textbf{Kanzara}}& \multicolumn{3}{c}{\textbf{Shirapur}}');
            elseif VILLAGE==1
                matrix2latex_newnew(Moments_4_table,[path1 'moments4_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{3}{c}{\textbf{Aurepalle}}');
            end
            % parameters chosen to target 4 moments and measurement error
            if pres==1
                rowLabels = {'n','$\delta$','s.e.','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\Delta V_{dc}$', '$\Delta \beta_{dc}$','$\frac{Var_{d\epsilon}}{Var_{dc}}$','s.e.','$\frac{Var_{d\nu}}{Var_{dy}}$','s.e.', 'Fit'};
            else
                rowLabels = {'n','$\delta$','s.e.','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\Delta V_{dc}$', '$\Delta \beta_{dc}$','$\frac{Var_{d\epsilon}}{Var_{dc}}$','s.e.','$\frac{Var_{d\nu}}{Var_{dy}}$','s.e.', 'Model criterion'};
            end
            if Nincproc>1
                matrix2latex_newnew(Moments_4_ms_table, [path1 'moments4ms_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{3}{c}{\textbf{Aurepalle}}& \multicolumn{3}{c}{\textbf{Kanzara}}& \multicolumn{3}{c}{\textbf{Shirapur}}');
            end
            % parameters chosen to target 4 moments and groupsize
        end
    end
    
    if Rov_sb==1 %|| one_per_aut==0
        columnLabels = {'Data','CD','Data','CD','Data','CD'};
        rowLabels = {'n','$\delta$','s.e.','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy\leq0}$', 'Model criterion'};
        matrix2latex_newnew(Moments_4_table,[path1 'moments4_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{2}{c}{\textbf{Aurepalle}}& \multicolumn{2}{c}{\textbf{Kanzara}}& \multicolumn{2}{c}{\textbf{Shirapur}}');
    end
    if one_per_aut==1 %|| one_per_aut==0
        columnLabels = {'Data','CD','Data','CD','Data','CD'};
        rowLabels = {'n','$\delta$','s.e.','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy\leq0}$', 'Model criterion'};
        matrix2latex_newnew(Moments_4_table,[path1 'moments4_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{2}{c}{\textbf{Aurepalle}}& \multicolumn{2}{c}{\textbf{Kanzara}}& \multicolumn{2}{c}{\textbf{Shirapur}}');
    end
    if sum(COALDEV-2==0)==1
        if isnan(phicalib)==1
            columnLabels = {'Data','SI','Data','SI','Data','SI'};
            rowLabels = {'n','$\delta$','s.e.','$\sigma$','\phi','s.e.','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy\leq0}$','$\frac{Var_{d\epsilon}}{Var_{dc}}$','s.e.','$\frac{Var_{d\nu}}{Var_{dy}}$','s.e.', 'Model criterion'};
            if Nincproc>1
                matrix2latex_newnew(Moments_4_ms_table_SI, [path1 'moments4ms_SI_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{2}{c}{\textbf{Aurepalle}}& \multicolumn{2}{c}{\textbf{Kanzara}}& \multicolumn{2}{c}{\textbf{Shirapur}}');
            end
            rowLabels = {'n','$\delta$','s.e.','$\sigma$','\phi','s.e.','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy\leq0}$', 'Model criterion'};
            matrix2latex_newnew(Moments_4_table_SI, [path1 'moments4_SI_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{2}{c}{\textbf{Aurepalle}}& \multicolumn{2}{c}{\textbf{Kanzara}}& \multicolumn{2}{c}{\textbf{Shirapur}}');
            matrix2latex_newnew(Moments_2_table_SI, [path1 'moments2_SI_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{2}{c}{\textbf{Aurepalle}}& \multicolumn{2}{c}{\textbf{Kanzara}}& \multicolumn{2}{c}{\textbf{Shirapur}}');
        end
        if COALDEV(1)==0
            columnLabels = {'Data','ID', 'CD','SI','Data','ID', 'CD','SI','Data','ID', 'CD','SI'};
        else
            columnLabels = {'Data', 'CD','ID','SI','Data','CD','ID','SI','Data','CD','ID','SI'};
        end
        if isnan(phicalib)==1
            if pres==0
                rowLabels = {'n','$\delta$','s.e.','$\sigma $','$\phi$','s.e.','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy\leq0}$', 'Model criterion'};
            else
                rowLabels = {'n','$\delta$','s.e.','$\sigma $','$\phi$','s.e.','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$','$\Delta V_{dc}$', '$\Delta \beta_{dc}$', 'Fit'};
            end
        else
            if pres==0
                rowLabels = {'n','$\delta$','s.e.','$\sigma $','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy\leq0}$', 'Model criterion'};
            else
                rowLabels = {'n','$\delta$','s.e.','$\sigma $','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$','$\Delta V_{dc}$', '$\Delta \beta_{dc}$', 'Fit'};
            end
        end
        if ext==1
            matrix2latex_newnew(Moments_4_table_plus, [path1 'moments4_plus_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{4}{c}{\textbf{Aurepalle}}& \multicolumn{4}{c}{\textbf{Kanzara}}& \multicolumn{4}{c}{\textbf{Shirapur}}');
            matrix2latex_newnew(Moments_2_table_plus, [path1 'moments2_plus_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{4}{c}{\textbf{Aurepalle}}& \multicolumn{4}{c}{\textbf{Kanzara}}& \multicolumn{4}{c}{\textbf{Shirapur}}');
            if pres==1
                rowLabels = {'n','$\delta$','s.e.','$\sigma $','$\phi$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\Delta V_{dc}$', '$\Delta \beta_{dc}$','$\frac{Var_{d\epsilon}}{Var_{dc}}$','s.e.','$\frac{Var_{d\nu}}{Var_{dy}}$','s.e.', 'Fit'};
            else
                rowLabels = {'n','$\delta$','s.e.','$\sigma $','$\phi$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\Delta V_{dc}$', '$\Delta \beta_{dc}$','$\frac{Var_{d\epsilon}}{Var_{dc}}$','s.e.','$\frac{Var_{d\nu}}{Var_{dy}}$','s.e.', 'Model criterion'};
            end
            if Nincproc>1
                matrix2latex_newnew(Moments_4_ms_table_plus, [path1 'moments4_ms_plus_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{4}{c}{\textbf{Aurepalle}}& \multicolumn{4}{c}{\textbf{Kanzara}}& \multicolumn{4}{c}{\textbf{Shirapur}}');
                columnLabels_plus = {'Data', 'CD','CD','ID','SI','Data','CD','CD','ID','SI','Data','CD','CD','ID','SI'};
                matrix2latex_newnew(Moments_4_ms_table_plus_plus, [path1 'moments4_ms_plus_plus_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels_plus, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{4+(countadd>0)}{c}{\textbf{Aurepalle}}& \multicolumn{4+(countadd>0)}{c}{\textbf{Kanzara}}& \multicolumn{4+(countadd>0)}{c}{\textbf{Shirapur}}');
           end
        end
    end
    if sum(COALDEV-3==0)==1
        columnLabels = {'Data','CD','Data','CD','Data','CD'};
        rowLabels = {'n','$\delta_{min}$','$\delta_{max}$','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy\leq0}$', 'Model criterion'};
        matrix2latex_newnew(Moments_4_table_FI, [path1 'moments4_FI_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{2}{c}{\textbf{Aurepalle}}& \multicolumn{2}{c}{\textbf{Kanzara}}& \multicolumn{2}{c}{\textbf{Shirapur}}');
    end
    
    if sum(COALDEV==0)==1 && sum(COALDEV-4==0)==1
        columnLabels = {'Data','ID', 'ID Het','Data','ID', 'ID Het','Data','ID', 'ID Het'};
        if pres==1
            rowLabels = {'n','$\delta$','s.e.','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\Delta V_{dc}$', '$\Delta \beta_{dc}$', 'Fit'};
        else
            rowLabels = {'n','$\delta$','s.e.','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy\leq0}$', 'Model criterion'};
        end
        matrix2latex_newnew(Moments_2_table_het,[path1 'moments2_' int2str(outlier) '_het' int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{3}{c}{\textbf{Aurepalle}}& \multicolumn{3}{c}{\textbf{Kanzara}}& \multicolumn{3}{c}{\textbf{Shirapur}}');
        % parameters chosen to target 4 moments
        matrix2latex_newnew(Moments_4_table_het,[path1 'moments4_' int2str(outlier) '_het' int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{3}{c}{\textbf{Aurepalle}}& \multicolumn{3}{c}{\textbf{Kanzara}}& \multicolumn{3}{c}{\textbf{Shirapur}}');
    end
    if est_rho1==1
        
        if COALDEV(1)==1
            columnLabels = {'Data', 'CD','ID','Data','CD','ID','Data','CD','ID'};
        else
            columnLabels = {'Data','ID','Data','ID', 'Data','ID'};
        end
        rowLabels = {'n','$\delta$','s.e.','$\sigma$','$\frac{Var_{dc}}{Var_{dy}}$','$\beta_{dcdy}$', '$\frac{Var_{dc|dy>0}-Var_{dc|dy\leq 0}}{Var_{dy}}$', '$\beta_{dcdy>0}-\beta_{dcdy\leq0}$','$\rho$','s.e.', 'Model criterion'};
        matrix2latex_newnew(Moments_4_rho1_table, [path1 'moments4rho1_cd0_' int2str(outlier) int2str(fixedeffect) Date '.txt'], 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', tablefontsize,'headline','& \multicolumn{3}{c}{\textbf{Aurepalle}}& \multicolumn{3}{c}{\textbf{Kanzara}}& \multicolumn{3}{c}{\textbf{Shirapur}}');
    end
end

